#include "PrimaryGeneratorAction.hh"

#include <G4Types.hh>

#include "DetectorConstruction.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleGun.hh"
#include "G4SystemOfUnits.hh"
#include "G4Threading.hh"
#include "MyMacros.hh"
#include "PrimaryGeneratorMessenger.hh"
#include "Randomize.hh"

// #include "G4ParticleTable.hh"
// #include "G4PhysicalConstants.hh"

constexpr G4double _beamEnergyDft{480.0};

PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* dtctr) : _fDetector(dtctr)
{
#if 0
    MYG4TID(tid);
    DINFOLN("```[tid {} PrimaryGeneratorAction/PrimaryGeneratorAction]", tid);
#endif

    InitializeMe();
}

void PrimaryGeneratorAction::InitializeMe()
{
    _fVerbose = _fDetector->GetVerbose();

    _fMessenger = new PrimaryGeneratorMessenger(this);
    _fParticleGun = new G4ParticleGun();

    _fPosition = G4ThreeVector(_fX0, _fY0, _fZ0);
    _fDirection = G4ThreeVector(0.0, 0.0, 1.0);

    SetBeamEnergy(_beamEnergyDft * MeV);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
    delete _fParticleGun;
    delete _fMessenger;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
    MYG4TID(tid);
#if 0
    DINFOLN("[tid {} PrimaryGeneratorAction] GeneratePrimaries()", tid);
#endif
    _fCounter++;

    // Simulation of beam position
    G4double x = _fX0;
    G4double y = _fY0;
    G4double z = _fZ0;  //2022.10.18
    _fPosition = G4ThreeVector(x, y, z);
    _fParticleGun->SetParticlePosition(_fPosition);

    do
    {  ///xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 单个点源
        if (_fDisTypeXY == "single")
        {
            // DINFOLN("[tid {} PrimaryGen/GeneratePrimaries] single case", tid)
            _fDirection = G4ThreeVector(0.0, 0.0, 1.0);
            _fDirection = _fDirection.unit();
            break;
        }

        ///xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx for QA beam model
        {
            G4double spoty = _fSpotDeltaY;
            // position of keypoint (between source and iso-center, in the spot-plane)
            G4ThreeVector posVSCenter = G4ThreeVector(0.0, spoty, 0.0);

            _fRMax2 = std::pow(5.0 * std::max(_fSigmaX, _fSigmaY), 2);  //2022.10.14
            G4double dx = 0, dy = 0;
            if (_fDisTypeXY == "flat")
            {
                dx = _fSigmaX * (1 - 2 * G4UniformRand());
                dy = _fSigmaY * (1 - 2 * G4UniformRand());
            }
            else if (_fDisTypeXY == "gauss")
            {
                do
                {
                    if (0.0 < _fSigmaX)
                    {
                        dx = G4RandGauss::shoot(0.0, _fSigmaX);
                    }
                    if (0.0 < _fSigmaY)
                    {
                        dy = G4RandGauss::shoot(0.0, _fSigmaY);
                    }
                } while (dx * dx + dy * dy > _fRMax2);
            }
            G4ThreeVector locPos = G4ThreeVector(dx, dy, 0.);
            G4ThreeVector dirX = G4ThreeVector(1.0, 0.0, 0.0);
            G4ThreeVector dirZ = (posVSCenter - _fPosition).unit();
            G4ThreeVector dirY = dirZ.cross(dirX);
            G4ThreeVector posVS = posVSCenter + RotationMatrix(locPos, dirX, dirY, dirZ);

            _fDirection = (posVS - _fPosition).unit();
            break;
        }
        ///xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx others

    } while (0);

    //DINFOLN("[tid {} PrimaryGen/GeneratePrimaries] momentum direction ({},{},{})", tid, _fDirection[0], _fDirection[1], _fDirection[2]);
    _fParticleGun->SetParticleMomentumDirection(_fDirection);

    ///000000000000000000000000000000000000000000000000 Simulation of beam kinetic energy
    G4double kinEnergy = _fEnergy;
    do
    {
        if (_fDisTypeE == "single")
        {
            kinEnergy = _fEnergy;
            break;
        }

        if (_fDisTypeE == "flat")
        {
            kinEnergy = _fEnergy + _fSigmaE * _fEnergy * (1 - 2 * G4UniformRand());
            break;
        }

        if (_fDisTypeE == "gauss")
        {
            if (0.0 < _fSigmaE) kinEnergy = _fEnergy + G4RandGauss::shoot(0.0, _fSigmaE * _fEnergy);

            break;
        }

    } while (0);

    _fParticleGun->SetParticleEnergy(kinEnergy);
    if (G4Threading::G4GetThreadId() == 0 && _fCounter == 1)
    {
        PrintSourceParameters();
    }

    if (_fVerbose > 0)
    {
        G4ParticleDefinition* particle = _fParticleGun->GetParticleDefinition();
        G4String particleName = particle->GetParticleName();
        LINFO(
            "```````````````````````````"
            "[PrimaryGeneratorAction/GeneratePrimaries] "
            "Event {}, Particle {}, PDGcode {}, KinEngy(MeV) {:>11.6F}, "
            "(cm) x {}, y {}, z {}",
            _fCounter, particleName, _fEnergy / MeV, x / cm, y / cm, z / cm);
    }

    _fParticleGun->GeneratePrimaryVertex(anEvent);
}

void PrimaryGeneratorAction::SetBeamEnergy(G4double val) { _fEnergy = val; }

void PrimaryGeneratorAction::PrintSourceParameters()
{
    MYG4TID(tid);
    G4ParticleDefinition* particle = _fParticleGun->GetParticleDefinition();
    G4String particleName = particle->GetParticleName();

    LINFO(
        "``````````````````````````````````````````````````"
        "[tid {} PrimaryGeneratorAction/PrintSourceParameters],\n"
        "Event {}, Particle {}, PDG code {},\n"
        "CentralEngy(MeV) {}, EngyDisType {}, sigma(rel.),\n"
        "CentralPosition(mm) ({};{};{}),\n"
        "PosDisType {}, Sigma(mm) ({};{};{}), RMax {},\n"
        "InitDirection ({};{};{}),\n"
        "CentralTheta {}, ThetaDisType {}, Sigma {}",
        tid, _fCounter, particleName, particle->GetPDGEncoding(),                           //
        _fEnergy / MeV, _fDisTypeE, _fSigmaE,                                               //
        _fX0 / mm, _fY0 / mm, _fZ0 / mm,                                                    //
        _fDisTypeXY, _fSigmaX / mm, _fSigmaY / mm, _fSigmaZ / mm, std::sqrt(_fRMax2) / mm,  //
        _fDirection.x(), _fDirection.y(), _fDirection.z(),                                  //
        0.0, _fDisTypeTheta, _fSigmaTheta / rad                                             //
    );
}

/// 逆变换矩阵，从 局域坐标系-> 全局坐标系
inline G4ThreeVector PrimaryGeneratorAction::RotationMatrix(G4ThreeVector& locPos,
                                                            G4ThreeVector& Rotx,
                                                            G4ThreeVector& Roty,
                                                            G4ThreeVector& Rotz) const
{  //2023.6.20
    G4double x = locPos.x() * Rotx.x() + locPos.y() * Roty.x() + locPos.z() * Rotz.x();
    G4double y = locPos.x() * Rotx.y() + locPos.y() * Roty.y() + locPos.z() * Rotz.y();
    G4double z = locPos.x() * Rotx.z() + locPos.y() * Roty.z() + locPos.z() * Rotz.z();
    return G4ThreeVector(x, y, z);
}

///0000000000000000000000000000000000000000000000000000000000000 backup
#if 0
  // Simulation of beam direction
  G4double ux = fDirection.x();
  G4double uy = fDirection.y();
  G4double uz = fDirection.z();

  // Beam particles are uniformly distributed over phi, cosTheta
  if(fDisTypeTheta == "flat") {
      G4double minCosTheta = std::cos(fSigmaTheta);
      uz = minCosTheta + (1.0 - minCosTheta)*G4UniformRand() ;
      ux = std::sqrt((1.0 - uz)*(1.0 + uz)) ;
  } else if(fDisTypeTheta == "gauss") {
      if(0.0 < fSigmaTheta){
          uz = std::cos(G4RandGauss::shoot(0.0,fSigmaTheta));
          ux = std::sqrt((1.0 - uz)*(1.0 + uz));
      }
  }

  G4double phi = twopi*G4UniformRand() ;
  uy = ux;
  ux *= std::cos(phi) ;
  uy *= std::sin(phi) ;
  fDirection = G4ThreeVector(ux,uy,uz) ;
#endif